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ABSTRACT 

The recent IceCube discovery of 0.1 - 1 PeV neutrinos of astrophysical origin opens up a new 
era for high-energy astrophysics. Although there are various astrophysical candidate sources, a 
firm association of the detected neutrinos with one (or more) of them is still lacking. A recent 
analysis of plausible astrophysical counterparts within the error circles of IceCube events 
showed that likely counterparts for nine of the IceCube neutrinos include mostly BL Lacs, 
among which Mrk 421. Motivated by this result and a previous independent analysis on the 
neutrino emission from Mrk 421, we test the BL Lac-neutrino connection in the context of 
a specific theoretical model for BL Lac emission. We model the spectral energy distribution 
(SED) of the BL Lacs selected as counterparts of the IceCube neutrinos using a one-zone 
leptohadronic model and mostly nearly simultaneous data. The neutrino flux for each BL Lac 
is self-consistently calculated, using photon and proton distributions specifically derived for 
every individual source. We find that the SEDs of the sample, although different in shape and 
flux, are all well fitted by the model using reasonable parameter values. Moreover, the model- 
predicted neutrino flux and energy for these sources are of the same order of magnitude as 
those of the IceCube neutrinos. In two cases, namely Mrk 421 and H 1914-194, we find a 
suggestively good agreement between the model prediction and the detected neutrino flux. 
Our predictions for all the BL Lacs of the sample are in the range to be confirmed or disputed 
by IceCube in the next few years of data sampling. 

Key words: astroparticle physics - neutrinos - radiation mechanisms: non-thermal - galax¬ 
ies: BL Lacertae objects: individual: Mrk 421, 1ES 1011+496, PG 1553+113, H 2356-309, 
1H 1914-194,1RXS J054357.3-553206 


1 INTRODUCTION 

During the first three years of data sampling in full configuration, 
the IceCube Neutrino Telescope based at the South Pole has reg¬ 
istered a sample of ver y high-energy neutrinos (3 0 TeV - 2 PeV) 
of astrophysical origin dAartsen et alJl2013L 120141 ). With a sample 
of 37 events the background only hypothesis with no astrophysical 
component has been rejected with a significance of more than 5 cr. 
A new era for high-energy astrophysics just started. 

Ex cluding, in this context, the possible co nnection with dark 
matter l lCherrv et alJl2014l : lEsmaili et alj|2014l) . the IceCube dis¬ 
covery calls for astrophysical sites, where cosmic-ray acceleration 
at the “knee” energy scale and above is efficiently at work, and 
the target, be it gas or photons, is sufficiently dense for interac- 
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tions with the cosmic rays. A firm association between the de¬ 
tected neutrinos and a specific class (or classes) of astrophysical 
sources is still lacking, notwithstanding the numerous scenarios 
that have been proposed so far. Models suggesting a pure Galactic 
origin of the observed neutrino signal can be strongly const rained 
by diffuse TeV-P eV y-ray limits (see e.g. lAhlers & Murasell2014l : 
Uoshi et alj |2014| and references therein), given that the produc¬ 
tion of very high energy (VHE) y-rays is an inevitable outcome 
of the hadronic interactions that lead to the neutrino production in 
the first place. However, a Galactic component contributing to the 
IceCube neutrino flux cannot be excluded at the moment. Proton- 
proton ( pp ) collisions in galaxy groups and/or star forming galax- 
ies have also been invoked to explain the diff use neutrino flux (e.g. 
[Loeb & Waxmard [20061 : M urase et aD 120131 ). Other extragalactic 
scenarios that include neutrino production through photohadronic 


(p7r) i nteractions in ac t ive galactic nuclei (AGN) (jS tecker et al 


1 199 it lMannheim|[l995l : lHalzen & Zaslll997l : iRachen & Meszaros 
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19981 

; Miicke et al. 19991; Atovan & Dermer 200 ll; Miicke et al. 

2003 

) and gamma-ray bursts (GRBs) ( 

Waxman & Bahcal 

1997; 

Murasej 

200i; Zhang & Kumar 2013 ; 

Asano & Meszaros 

2014; 

Revnoso 

20141 Petronoulou et alj 2014 

Winter et al, 20141) have 


been extensively discussed in the literature. 

A commonly adopted approach for the calculation of the neu¬ 
trino flux at Earth (see e.g. Munnhejmetal.| [200 ll : iBaerwald et alj 
120 1 3l : iKTstler et alJl2013l : lMurase et all20l4) is the following: (i) a 

generic distribution of accelerated protons in the source is assumed; 
(ii) a generic photon distribution that acts as a target field for pn 
interactions is also assumed; (iii) the neutrino spectrum is then cal¬ 
culated; and (iv) the neutrino flux is integrated over redshift assum¬ 
ing a luminosity function for the particular astrophysical source. If 
many sources of the same type contribute to the observed flux, then 
the approach described above is the most relevant one. A differ¬ 
ent approach is to focus on a few individual sources of the same 
class and derive the respective neutrino fluxes using proton and 
photon distributions that are un ique to each object dMiicke et alj 
120031 ; iDimitrakoudis et al.ll20l4) . In this regard, information about 
the directions of the neutrino events on the sky is important. This 
approach may still lead to a study of the diffuse neutrino flux, from 
a certain class of sources, but with the caveat that it may not repre¬ 
sent the total diffuse flux; the latter could be arrived at only upon an 
estima te of the percentage of that class’s contribution to the total. 

IPadovani & Resconil d2014l ) - henceforth PR 14, have recently 
searched for plausible astrophysical counterparts within the error 
circles of IceCube neutrinos using high-energy y -ray (GeV-TeV) 
catalogs. Assuming that each neutrino event is associated with one 
astrophysical source and using a model-independent method they 
derived the most probable counterparts for 9 out of the 18 neutrino 
events of their sample. Interestingly, these include 8 BL Lac ob¬ 
jects, amongst which the nearest blazar, Mrk 421, and two pulsar 
wind nebulae. Although Flat Spectrum Radio Quasars (FSRQs) are 
belie ved to be more effic ient PeV neutri no emitters than BL La cs 
(e.g. iDermer et aid 120141 ; see, however, iTavecchio et aid d 201 4| ) ). 
these do not appear in the list of most probable counterparts. This 
is not necessarily related to the physics of neutrino emission from 
different types of blazars but could be a consequence of the “ener¬ 
getic” test applied by PR14 (for more details, see §4 therein). 

Motivated by the aforementioned results, we investigate the 
BL Lac-PeV neutrino connection within a specific theoretical 
framework for blazar emission where the y -rays are of photo- 
hadronic origin. In contrast to the majority of works on neutrino 
emission from blazars, we do not assume a generic SED but we fit 
instead the multiwavelength (MW) photon spectra using the radi¬ 
ation of particles whose distributions have been accordingly mod¬ 
ified by the cooling and escape processes. In our model the low- 
energy emission of the blazar SED is attributed to synchrotron ra¬ 
diation of relativistic electrons, whereas the observed high-energy 
(GeV-TeV) emission is the result of synchrotron radiation from 
pairs produced by charged pion decays. Pions are the by-product 
of p7t interactions of co-accelerated protons with the internally pro¬ 
duced synchrotron photons. Our approach allows us to identify the 
observed y-ray emission as the emission from the p n component 
and associate it with a very high-energy (~ 2-20 PeV) neu¬ 
trino emission. Analytical expressions for the typical energies of 
neutrinos and photons produced through photohadronic processes, 
the proton cooling rate due to photopion and photopair interac¬ 
tions, as well as the dependence of all aforementioned quanti- 
ties on o bservable q uantities of blazar emission can be found in 
IPetropoulou & MastichiadiU j2015i) . 

First results of the model were presented in 


iDimitrakoudis et all J20 1 4l) - henceforth DPMI4, where the 
neutrino flux expected from the prototype blazar Mrk 421 in a low 
state was derived. The prediction of the model muon neutrino flux 
was tantalizing l v clos e to the IC-40 sensitivity limit calculated by 
iTchernin et all 120131) for Mrk 421 (see Fig. 1 in DPM14). We 
extend our previous analysis to include those BL Lacs identified as 
the most probable counterparts of the IceCube neutrinos by PR 14 
(see Table 4 therein). Our aim is threefold: (i) to examine whether 
the working model (the so-called LH/r model, see DPM14) can fit 
the SED of the aforementioned blazars - there was no guarantee 
before our attempts, since we did not select a priori sources with 
specific spectral characteristics; (ii) to calculate the neutrino flux 
relative to the y-ray one, should acceptable fits be obtained, and 
(iii) to compare the model-derived neutrino fluxes with those 
calculated for each single neutrino event. We show that Mrk 421 
does not constitute a unique case. Instead, all BL Lacs in our 
sample can be modelled by the leptohadronic model and for most 
of them the model-predicted neutrino flux is close to the observed 
one, at least within the 3<x error bars. 

This paper is structured as follows. We start with a short de¬ 
scription of the theoretical model and the fitting method we used 
in Sect. [2] In Sect. [3] we provide general information about the ob¬ 
servational datasets used for each BL Lac. In Sect. [4] we present 
the combined MW photon and neutrino spectra for each of the case 
studies, focusing on the ratio of the model-predicted neutrino lu¬ 
minosity to the y-ray luminosity and its dependence on observable 
quantities. We also comment on the model-derived energetics for 
each source. We discuss our results in Sect.[5]and conclude with a 
summary in Sect. [6] Throughout this work we have adopted a cos¬ 
mology with Q m = 0.3, Da = 0.7 and H 0 = 70 km s' 1 Mpc~*. 
For the attenuation of the VHE part of the spectra due to photon- 
photon pair productio n on the extragalact i c back ground light (EBL) 
we used the model of lFranceschini et all d2008l) . 


2 MODEL AND METHOD 
2.1 Model description 

T he physical model we use h as been described, in general terms, 
in IDimitrakoudis et all d2012l) - henceforth D MPR12, and, in the 
contex t of modelling Mrk 421 in particular, in iMastichiadis et~al] 
d2013h and DPM14. We consider a spherical blob of radius R, con¬ 
taining a tangled magnetic field of strength B and moving towards 
us with a Doppler factor d. Protons and electrons are accelerated 
by some mechanism whose details lie outside the immediate scope 
of this work, and are subsequently injected isotropically in the vol¬ 
ume of the blob with a constant rate. Their interaction with the 
magnetic field and with secondary particles leads to the develop¬ 
ment of a system where five stable particle populations are at work: 
protons, which lose energy by synchrotron radiation, Bethe-Heitler 
(pe) pair production, and photopion (pn) interactions; electrons (in¬ 
cluding positrons), which lose energy by synchrotron radiation and 
inverse Compton scattering; photons, which gain and lose energy 
in a variety of ways; neutrons, which can escape almost unimpeded 
from the source region, with a small probability of photopion inter¬ 
actions; and neutrinos, which escape completely unimpeded. 

The interplay of the processes governing the evolution of 
the energy distributions of those five populations is formulated 
with a set of time-dependent kinetic equations. Through them, en¬ 
ergy is conserved in a self-consistent manner, since all the en¬ 
ergy gained by a particle type has to come from an equal amount 
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of energy lost by another particle type. To simultaneously solve 
the coupled kinetic equations for all particle types we use the 
time-dependent code described in DMPR12. Photopion interac¬ 
tions are modelled using the results of the Monte Carlo event 
generator SOPHIA jMucke et alJbood) . while Bethe-Heitler pair 
production is similar l y mod elled with the Monte Ca r lo resu lts of 
IProtheroe & Johnsonl ll 1 9961) and iMastichiadis et all 120051) . The 
only particles not modelled with kinetic equations are muons, 
p ions and kaons. Their treatment is described in DPMI4 and 
IPetropoulou et alj d2014l) but is of negligible importance for the 
range of parameter values used in the present work. 

The prime mover of all processes is the rate at which energy 
is injected in relativistic protons and electrons. This, in terms of 
compactnesses, is expressed as: 


r .. . - 


L' .<x T 

i,inj 1 

AnRm\ri ’ 


0 ) 


where i denotes protons or electrons (i=e,p), cr T is the Thomson 
cross section, and L' jnj denotes the injection luminosity as mea¬ 
sured in the rest frame of the emitting region. In principle, the ra¬ 
dius R and Doppler factor S can be constrained using variability 
arguments. However, we prefer to treat them as free parameters, 
since information for high-energy variability is not available for all 
sources in the sample. The escape times for protons and electrons 
are assumed to be the same (t p>esc = t ee sc), and complete the set 
of free parameters dependent on the geometry and kinematics of 
the blob. What remains to be defined are the shapes of the proton 
and primary electron injection functions, which are described as 
power-laws with index jj extending from yj im in to yi, max . In some 
cases, modelling of the SED requires the primary electron distri¬ 
bution to be a broken power-law, with the break energy y ebr and 
the power-law index above the break ,s Ci h being two additional free 
parameters. 

As long as the escaping protons and neutrons are energetic 
enough, they are susceptible to photopion interactions with ambi¬ 
ent photons in Galactic and intergalactic space, such as the cos¬ 
mic microwave and infrared b ackgro unds, producing additional 
high-energy neutrinos (iStecken [l973|) . Neutrons also rapidly de¬ 
cay into protons I Sikora et all 19871: Kirk & Mastichiadlsl 1 19891 : 


iGiovanoni & Kazanasll99d : lAtovan & Dermeil200 ll) . leading also 

to neutrino production. In this work we focus on the neutrino emis¬ 
sion from the objects themselves, and will not consider any such 
additional contributions from escaping high-energy nucleons. 

The reader should note that we placed our emphasis on fitting 
the MW data with a variant of the leptohadronic model where the 
y-ray emission is attributed to synchrotron radiation of the secon¬ 
daries arising from charged pion decay. Clearly, this is not the only 
way of fitting the SED, as it is possible to fi t the y-rays with proton 
synchrotron radiation (e.g. lAharoniarfeOOOl) . In that case, however, 
the model produces a very low neutrino flux at ~EeV energies (see 
LHs model in DPM14), which makes it irrelevant to the IceCube 
TeV-PeV detections. Moreover, scenarios where the high-energy 
blazar emission is the result of y-rays produced through pp col¬ 
lisions within the blazar jet, constitute another alternativqj. Such 
models would require, however, uncomfortably high densities of 
non -relativistic (cold) proton s to be present in the blazar jet (see 
e.g. lAtovan & Dermeifoool and references therein). Last but not 
least, there is always the possibility of assuming a purely leptonic 


1 To our knowledge, direct application of such models to MW observations 
of blazars is still lacking in the literature. 
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Figure 1. Flow diagram illustrating the steps of the fitting method. 


model feottcher et al J20 13l~) which, however, does not produce any 
neutrinos. 


2.2 Fitting method 

The multi-parameter nature of the problem (11-13 free parameters; 
see Sect. 0 and the long simulation time, which is required for 
every parameter set (of the order of ~ 0.5 d), prevents us from 
a blind search of the whole available parameter space. Moreover, 
the DMPR12 code in its current form has not been developed for 
parallel computing. Yet, the application of the model to actual ob- 
servations is still viable (e.g.I Mastic hia dis et aLl l2013l ; |Petropouloul 
l2014h thanks to the insight we gain from an analytical treatment of 
the problem. Thus, instead of automatically scanning the parameter 
space, we start the fitting procedure with a set of values, which are 
physically motivated (see below). Then, we perform a series of nu¬ 
merical simulations with parameter values lying close to the initial 
parameter set, until a reasonably good fit to the SED is obtained. 
The term “reasonably good” fit refers to a model SED that passes 
through most observational points. This is not necessarily the best 
fit in the strict sense, i.e. the one with the lowest^ 2 value. However, 
the main results of our study, such as the y-ray/neutrino connection, 
are not affected by the exact value of the^ 2 statistics. 

The steps of the fitting method are illustrated in Fig. |T| The 
parameters of the problem can be divided into two groups; those 
related to the properties of the emission region, i.e. B, R and S and 
those related to the primary particle distributions (see two columns 
in Fig. [T}. As benchmark values for the magnetic field we chose 
B = 0.1,1,10 G. These cover the range of low-to-moderate mag¬ 
netic fields. Higher values, e.g. 6 > 20 G, were not included to the 
benchmark set, since they tend to favour proton synchrotron radia- 


els, see Aharonian 2000j; Miicke et al. 20031; 

Bottcher et al. 20131; 

Dimitrakoudis et al. 20141; Cerruti et al. 2014 

). For each value of 


B, we then adopted R = 10 15 cm and 10 16 cm and the set of bench¬ 
mark values is completed with (5=15 and 30. We did not search in 
the first place for fits having even higher values of 6 and R, since 
they would decrease the optical depth ( f w ) for photopion interac¬ 
tions and for neutrino production. If the low-energy hump of the 
SED has peak spectral index ft, peak frequency v s and total lumi¬ 
nosity L syn the optical depth for neutrino production is given by 
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(e.g. PM15) 


/p*(y p ) “ 4.4 x 10- 


4yn,45 4(/L V s ) 


^15<5j V s> 16 ( 1 + z) 


2 7 p 


Jri 
V/p 


,7p>rr 


The characteristic photon energy is given by 
B , 


2v,„ 


( 2 ) 


where we introduced the notation q x = ry/ l 0 X in cgs units, A is 
a function of ji and v s (see equations (23) and (27) in PM15), z 
is the redshift of the source, v min is the low cutoff frequency of the 
low-energy spectrum (for a continuous power-law), y p is the proton 
Lorentz factor and y p p,r) is the minimum Lorentz factor that satisfies 
the threshold condition for pn interactions with the peak energy 
photons. This is given by (for details, see PM15) 


D C r 


where 

7e,p^ 


V 4 m e rp ’ V_ °- 2 


( 10 ) 


( 11 ) 


Using equations 0-0. Q- and 0 we derive a rough estimate of 
the proton injection compactness: 

3 f-y.46 + z) 


fin ini — 2.3 X 10 

PJ L syn ,45<5i VMfi, V.) 


( 12 ) 


7 p plr> ~ 3.5 x 10 7 <5iv-J 6 (l + z) -1 . 


(3) 


Protons with the aforementioned Lorentz factor lead to the produc¬ 
tion of neutrinos with energy, as measured in the observer’s frame, 
given by 


e v =4 17.5 PeV<5|(l +z) _2 v" 


(4) 


/ V S ,16(1 + Z)\ 1A 


need a stronger magnetic field (for modelling of blazar flares, see 

(5) 

Mastichiadis & Kirk 1997^Coppi & Aharonian 1999i;[Sikora et alJ 

\ <5iB-i / 


200ll;lKrawczvnski et alj 20o3). which in turn would affect the re- 


In most cases, the above expressions is a safe estimate of the peak 
energy of the neutrino spectrum. 

We use equation 0 as the minimum requirement for the max¬ 
imum proton energy of the proton distribution, i.e. y Pima x = A-y p p,r) 
with k a numerical factor of ~ 2-5. On the other hand, the maxi¬ 
mum (or break) Lorentz factor of the primary electron distribution 
is related to the peak frequency of the low-energy hump as 


7e,max /br — 6 X 10 


For the lower cutoffs of the proton and primary electron injection 
functions we use as default values 7 e ,pmin = 1 (grey colored param¬ 
eters in Fig.0. If we find large deviations of the model SED from 
the observations at v < v s , we search for y emin > L We search for 
different y ppnin than the default one only if the proton distribution 
has s p > 2 and extension to y p , m i n = 1 leads to extremely high 
injection luminosities. We note, however, that this is not a typical 
case. 

The power-law index of the primary electron distribution (i e 
or i e h ) is directly related to the observed spectral index of the low- 
energy component in the SED. However, it is not straightforward 
to relate j p with the observed photon index, e.g. in the Fermi/LAT 
regime, where the emission from secondaries produced in photo- 
hadronic interactions dominates. Thus, as a first estimate we use 

5 p — Sq. 

The injection luminosity of primary electrons is also directly 
related to the observed lum inosity of the low-energy hump as (e.g. 
iMastichiadis & Kirl J 19971) 


Summarizing, for each set of benchmark values for B,R and S and 
equipped with the analytical estimates of equations Ol d 12b and 
observables, we may significantly reduce the parameter space we 
need to scan (Fig.0. 


3 CASE STUDIES 

Our modelling targets are BL Lacs that are spatially and energet¬ 
ically correlated with neutrino detections. BL Lacs are known for 
their variable emission at almost all energies. In high-states they 
usually exhibit major flares, especially in X-rays and y-rays. Mod¬ 
elling of their broadband emission in low and high states requires 
a different parameterization. For example, a high-flux state might 


4,inj = 2 x 10 


-4 ^syn,45 


( 6 ) 


To derive a first estimate of the proton injection luminosity we as¬ 
sume that the y-ray (> 10 GeV) luminosity is totally explained by 
the synchrotron emission of p;r pairs: 

l 7 « (7) 

W<) - 

where f prr is given by equation 0. The factor 1/8 in equation 0 
comes from the assumption that half of the produced pions produce 
electrons/positron pairs, and each of them gains ~ 1/4 of the parent 
proton energy: 

„2 


e p = y p m p c- 


(9) 


suiting y-ray and neutrino fluxes. Thus, it is important to build an 
SED using simultaneous data or, if this is not possible, using MW 
observations that have been obtained at least within the same year. 

Out of the 8 BL Lacs reported by PR 14, two have unknown 
redshifts: SUMSS J014347-584550 and PMN J0816-1311, which 
are the probable counterparts of neutrinos 20 and 27, respectively. 
Since redshift is a crucial parameter in determining the luminos¬ 
ity of the emitting regions and, thus, the interplay of processes at 
work there, we have chosen to exclude those two BL Lacs from this 
study. Therefore, we are left with six targets, which are: 

• Mrk 421 at z = 0.031 dSbarufatti et al.|[2005t) . This is the clos¬ 
est high-frequency peaked blazar (HBL) and the target of many 
MW campaigns. Mrk 421 was also the first source which our model 
(see Sect. l2.lt was applied to and the main focus of an earlier work 
(DPM14). We first quote the results of DPM14, that used the simul¬ 
taneous X-ray (RXTE) and TeV y-ray (Whippl e and HEGRA) data 
of the 2001 MW campaign that lasted 6 days dFossati et alj|2008h . 
In particular, DPM14 fitted the SED of March 22nd/23rd 2001 that 
corresponds to a pre-flare state. During this period no GeV observa¬ 
tions were available. Th us, th e Fermi/LAT data during the period of 
the IC-40 configuration ( lAbdo etalj|201lh were included for com¬ 
parison reasons_ onlwjAs a second step, we use the MW observa¬ 
tions of lAbdo et alj d201 ll) together with the neutrino event ID 9 in 
order to build a hybrid SED. To avoid repetition and to emphasize 
the prediction of DPMI4 regarding the muon neutrino flux, we do 
not attempt to fit the new dataset. Instead, we use the same param¬ 
eters as in DPMI4, apart from the Doppler factor, and compare the 
model results against the [Abdo et al.1 d201 ill observations and the 
neutrino flux for the event ID 9. 

• PG 1553+113 at z = 0.4 jAleksic et al.l2012l) . This is not only 
the most distant source of this sample, but also belongs to the class 
of most distant TeV blazars in general. Its redshift has not been 
firmly measured, yet the value adopted here can be considered as a 
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safe e stimate (for different techniques and results, see lAleksic et alj 
|20I3) . Similarly to 1ES 1011+496 below, it has been only re¬ 
cently detected in V HE y-rays by both H .E.S.S. dAharonian et alj 
l2006bl) and MAGIC dAlhert et alj|2007^1 . The data from the sec¬ 
ond year Fermi-LAT Sources Catalog at ~ 1 GeV (2FGLJc) 0 
indicate a variable high-energy flux. The 2005-2009 MAGIC ob¬ 
servations also suggest varia bi lity at the VHE part of the spectrum 
(see Fig. 1 in lAleksic et alj|2012h . In 2005 the source exhibited 
large flux variations in X-rays, yet there were no simultaneous data, 
neither at GeV n or at Te V ener gies. We adopted the data shown 
in Fig. 6 of lAleksic et alj d2012h . For the SED fitting we used, in 
particular, the X-ray data that correspond to the intermediate flux 
level of 2005 along with the time-averaged 2008-2009 Fermi/LAT 
data and the MAGIC observations for the periods 2005-2006 and 
2007-2009. We include in our figures the high state X-ray data of 
2005 only for illustrative reasons. In principle, one could try to 
fit the time-averaged y-ray data together with the 2005 high state 
in X-rays. However, we find less plausible a scenario where the 
2005 X-ray high-state is not accompanied by an increased y-ray 
flux but, instead, is related to the average y-ray flux level. There 
are two reasons for this: the y-ray and X-ray e missio n in blazars 
are usually correlated (e.g. for Mrk 421, see iFossati et al.l l2008h 
and the y-ray emission of PG 1553+113 is itself variable (see e.g. 
lAharonian et all2006bl l. 

• 1ES 1011+496 at z = 0.212 dAlbert et alj|2007bh . This is an 

HBL source at indermediate redshift, that has been detected for 
the fir st time in TeV y-rays in 2007 with MAGIC l lAlbert et alj 
l2007bl) . The discovery was the result of a follow up observation 
of the optical high state reported by the Tuorla blazar monitor¬ 
ing prograrrE The previous non-detections in TeV y-rays suggest 
that 1ES 1011+496 is variable both in optical and VHE y-rays. 
Its X-ray emission i s described by a stee p spectrum in both low 
and high states (see Rc inthal et a l .1 1 20121 and references therein). 
Because of the variability observed across the MW spectrum, we 
chose to apply our model to simultaneous data, if possible. We 
therefore used data fro m the f irst MW campaign in 2008 as pre¬ 
sented in lReinthal et alj 12012l l. despite the fact that they are char- 
acteri z ed as prelimin ary. We note that the MAGIC data presented in 
Rcinthal et al J d2012l) ha ve been de-absorbed using the EBL model 
of Kneiske & Dolel <l201Ch . Using the same model we retrieved the 
observed MAGIC spectrum, which we show in the respective fig¬ 
ure. We also took into acc ount the Fermi/LAT data from the firs t 
(1FGL) l lAbdo et alj|2010l) and second (2FGL) dNolan et alj|2012l) 
catalogs. In the respective figures we include 2FGLJc data, in order 
to show the variability in the GeV regime. _ 

• H 2356-309 at z = 0.165 dFalomd fl99lh . This is one of 
the brightest HBLs in X -rays, whic h led to its early X-ray de¬ 
tection with UHURU flForman et alJI 19780 . Its X-ray emission is 
variable bo th in flux and in spectral shape (for more information, 
see E E.S.S. Collaboration et alJl20ldt . TeV y-ray emission from 
H 2356-309 was observed for the first time in 2004 by H.E.S.S. 
dAharonian et al Jl2006ah . Using the SED builder too@ of the ASI 
Science Data Centre (ASDC) we chose those X-ray observations 
that fall within the period 2004-2007 of H.E.S.S. observations. We 
also included data from the 1FGL and 2FGL catalogs, although 
they refer to the period 2008-2010, as well as hard X-ray data ob¬ 
tained by Swift/BAT between 2004 and 2010. Finally, we over- 


Dec. [deg] 



Figure 2. Sky map in equatorial coordinates of the five IceCube neutrino 
events (crosses) that have as probable astrophysical counterparts the six 
BL Lac sources (stars) mentioned in text. The red circles correspond to 
the median angular error (in degrees) for each neutrino event (see TableQ}. 


plotted the historical X-ray high state observed with BeppoSaX, 
for comparison reasons. __ 

• 1H 1914-194 at z = 0.137 JCarangelo et alj|2003lt . This is the 

first BL Lac of the sample that has not yet been detected in TeV 
y-rays, and the least well covered in different energy bands among 
the sources of this sample. In order to apply our model we con¬ 
structed a SED using available data from the ASDC, which is com¬ 
prised mainly of NED archival data. Given the quality of the avail¬ 
able dataset, our main focus was to describe as well as possible the 
Fermi/LAT data, while roughly describing the observations in other 
energy bands. _ 

* 1 RXS J054357.3-553206 at z = 0.374 istadnik & Romani! 
l2014h . This is the second most distant object of the sample, after 
blazar PG 1 553+ 11_3. It w as first detected and classified as a BL Lac 
by ROSA T dFischer et al.|| 19981 ). Following radio observations by 
lAnderson & Filipovic d2009h it was classified as an HBL, and it 
was s oon thereafter detected in y-rays by Fermi/LAT dAbdo et alj 
120101) . It rem ains still undetected in TeV energies. Its redshift was 
estimated by IStadnik & Romani J2014lf . based on the assumption 
that the flux of the host elliptical galaxies o f BL Lacs can be used 
as a standard candle, and also by lPita et al .1 ( 120141) . based on spec¬ 
troscopy on a weak Ca I I H/K doublet and a Na I absorption line. 

: Pita et al. ( 2014ll estimate (z = 0.237) is somewhat lower 


While the 


than the IStadnik & Roman] 120141) one, adopting it would not sig¬ 
nificantly alter out fit. Using the SED builder tool of ASDC we 
combined infrared (WISE), X-ray (Swift/XRT and XMM-Newton) 
and y-ray (Fermi/LAT) observations that were obtained within the 
same year (2010). We also included older XMM-Newton observa¬ 
tions performed in 2005 and 2007, for comparison reasons. Finally, 
the 2FGLJc data illustrate the variability in the GeV energy band. 


Before closing this section and for completeness reasons, we sum¬ 
marize in TableQ] the coordinates of the five IceCube neutrinos and 
the six most probable blazar counterparts, as well as their angular 
offset. We also included the median angular error for each neutrino 
detection and the respective detection time in Modified Julian Days 
(MJDfi A sky map (in equatorial coordinates) of the five IceCube 
neutrinos and their respective counterparts is shown in Figure[2] 


3 http://www.asdc.asi.it/fermi2fgl/ 

3 http://users.utu.fi/kani/lm 5 We note that all data listed in Table[]]are adopted from Tables 1 and 2 in 

4 http://www.asdc.asi.it/SED PR14. 
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Table 1. List of the most probable counterparts of selected high-energy IceCube neutrinos used in our modelling. 


IceCube ID 

Ra(deg) 

Dec(deg) 

Median angular error (deg) 

Time (MJD) 

Counterpart 

Ra(deg) 

Dec(deg) 

Angular offset (deg) 

9 

151.25 

33.6 

16.5 

55685.6629638 

Mrk 421 

166.08 

38.2 

12.8 






1ES 1011+496 

155.77 

49.4 

15.9 

10 

5.00 

-29.4 

8.1 

55695.2730442 

H 2356-309 

359.78 

-30.6 

4.7 

17 

247.40 

14.5 

11.6 

55800.3755444 

PG 1553+113 

238.93 

11.2 

8.9 

19 

76.90 

-59.7 

9.7 

55925.7958570 

1RXS J054357.3-553206 

85.98 

-55.5 

6.4 

22 

293.70 

-22.1 

12.1 

55941.9757760 

1H 1914-194 

289.44 

-19.4 

4.8 


all processes - v„ - Fermi 2011 i—e—i 

rtn- optical Whipple + HEGRA i—•—i 

v e +v^ . RXTE ■ IC-40 limit - 

logs (eV) 

-4 -2 0 2 4 6 8 10 12 14 16 


! 


I 
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Figure 3. Left panel: SED of blazar Mrk 421 as modelled in DPM14. The Fermi/LAT points are shown only for illustrative reasons, as they were not included 
in the fit of the original paper DPM14. Right panel: SED of Mrk 421 averaged over the period 19 January 2009 - 1 June 2009 (Abdo et al. 2011); all data 
points are from Fig. 4 in Abdo et al. (2011). The model SED (black line) and neutrino spectra (red lines) are obtained for the same parameters as in DPM14, 
except for the Doppler factor, i.e. 6 = 20. The neutrino event ID 9 is also shown. In both panels, the model spectra are not corrected for absorption on the EBL. 
For comparison reasons, the unattenuated y-ray emission from the n° decay is over-plotted with dash-dotted grey lines. 


4 RESULTS 

We first present the results of our SED modelling along with 
the predicted neutrino flux for each source, i.e. we build model 
hybrid SEDs. The parameter values derived for each source are 
summarized in Table [2] Various aspects of the adopted theo- 
retical framewor k were investigated in detail in a recent work 
dPetropoulou & Mastichiadisll2015l - hereafter PM15). For details 
and analytical expressions that give insight to the results that fol¬ 
low, we refer the reader to PM15. In the last two subsections we 
discuss in more detail the neutrino emission and the energetics for 
each source. 

4.1 SED modelling 

4.1.1 Mrk 421 

As already pointed out in the introduction, the most intriguing re¬ 
sult of the leptohadronic modelling of Mrk 421 was the predicted 
muon neutrino flux, which was close to the IC-40 sens itivit y limit 
calculated for the particular source bv lTchernin et al .1 (120131) . This 
is shown in Fig. [3] (left panel). We note that in the original pa¬ 
per DPM14,_whose results we quote here, the Fermi/LAT data 
dAhdo et ai]|201 ll) shown as grey open symbols in Fig. [3] were not 
included in the fit, since they were not simultaneous with the rest 
of the observations. The total and muon neutrino spectra are plot¬ 
ted with dotted and solid red lines, respectively. We use the same 
convention in all figures that follow, unless stated otherwise. The 
bump that appears two orders of magnitude below the peak of the 


neutrino spectrum in energy originates from neutron decay; we re¬ 
fer the reader to DPM14 for more details. 

According to PR 14, Mrk 421 is the most probable counterpart 
of the corresponding IceCube event (ID 9). This is now included 
in the right panel of Fig. [3] where the time-ave raged S ED during 
the MW campaign of January 2009 - June 2009 dAbdo et al.ll201 11) 
is also showitj. As we have already pointed out in Sect. [3] in or- 
d er to avo id repeti tion, we did not attempt to fit the new data set of 
lAbdo et all d201 ll) . Thus, the model photon (black line) and neu¬ 
trino (red line) spectra depicted in the right panel of Fig. [3] were 
obtained for the same parameters used in DPM14 (see also Tabled, 
except for the Doppler factor, which was here chosen to be 6 = 20. 
The fit to the 2009 SED is surprisinly good, if we consider the fact 
that we used the same parameter set as in DPM14 with only a small 
modification in the value of the Doppler factor. 

We note that the MW photon spectra shown in both panels 
have not been corrected for EBL absorption in order to demon¬ 
strate the relative peak energy and flux of the neutrino and n° com¬ 
ponents. The y-ray emission from the n° decay is shown in both 
panels as a bump in the photon spectrum, peaking at ~ 10’° Hz. Al¬ 
though the model prediction for the ratio of the n° y-ray luminosity 
to the total neutrino luminosity is roughly two (see e.g. DMPR12), 
the 7r° component shown in Fig. [3] is suppressed because of inter¬ 
nal photon-photon absorption. To demonstrate better the effect of 
the latter, we over-plotted the n° component that is obtained when 
photon-photon absorption is omitted (dash-dotted grey lines). In 

6 All data points are taken from Fig. 4 in Abdo et all J20 1 ll ). 
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Table 2. Parameter values for the SED modelling of blazars shown in Table 4 in PR 14. The ordering of presentation is the same as in Sect.0 The parenthesis 
below the name of each source encloses the ID of the associated neutrino event. The double horizontal line separates parameters used as an input to the 
numerical code (upper table) from those that are derived from it (lower table). 


Parameter 

Mrk 421 
(ID 9) 

PG 1553+113 
(ID 17) 

1ES 1011+496 
(ID 9) 

H 2356-309 
(ID 10) 

1H 1914-194 

(ID 22) 

1RXS J054357.3-553206 
(ID 19) 

z 

0.031 

0.4 

0.212 

0.165 

0.137 

0.374 

B (G) 

5 

0.05 

0.1 

5 

5 

0.1 

R (cm) 

3 x 10 15 

2 x 10 17 

3 x 10 16 

3 x 10 15 

3 x 10 15 

3 x 10 16 

6 

26.5 

30 

33 

30 

18 

31 

ye,min 

1 

1 

1 

1.2 x 10 3 

1.2 x 10 

1 

Te,max 

Bx 10 4 

2 x 10 6 

10 5a 

10 5 

10 5 

1.2 x 10 5a 

Te, br 

- 

2.5 x 10 5 

- 

2 x 10 4 

10 3 

- 

Sc, 1 

1.2 

1.7 

2.0 

1.7 

2.0 

1.7 

■?e,h 

- 

3.7 

- 

2.0 

3.0 

- 

f ■ ■ 
t'e.inj 

3.2 x 10- 5 

2.5 x 10- 5 

4 x 10- 5 

8 x 10~ 6 

6 x 10^ 5 

2.5 x 10- 5 

yp,min 

1 

10 3 

1 

1 

1 

1 

Tp,max 

3.2 x 10 6 

6 x 10 6 

1.2 x 10 7 

10 7 

2 x 10 6 

1.2 x 10 7 

5 P 

1.3 

2.0 

2.0 

2.0 

1.7 

2.0 

f ■ • 

^p/nj 

2 x 10~ 3 

3 x 10^ 4 

4 x 10~ 3 

2.5 x 10- 3 

io~ 2 

8 x 10~ 4 

L e,mj ( er g/ S > b 

4.4 x 10 40 

2.3 x 10 42 

5.6 x 10 41 

1.1 x 10 40 

8.3 x IO 40 

3.5 x 10 41 

L p,inj ( er g/ S > b 

5.1 x 10 45 

5.1 x 10 46 

10 47 

6.4 x 10 45 

2.5 x 10 46 

2 x 10 46 

L y .TeV (erg/s) c 

1.7 x 10 45 

7.9 x 10 46 

5 x 10 45 

5 x 10 44 

10 45 

6.3 x 10 45 

L v (erg/s) d 

2.5 x 10 45 

8.1 x 10 45 

2.5 x 10 45 

4 x 10 44 

2 x 10 45 

6.3 x IO 44 

Tv/ 

1.5 

0.1 

0.5 f 

0.8 

2.0 

0.1 

f P n B 

2 x 10’ 5 

6 x 10~ 6 

1.3 x 10~ 6 

3.1 x 10~ 6 

1.2 x 10~ 5 

1.1 x 10- 5 


a An exponential cutoff of the form e~ y l yt ' max was used. 
b Proton and election injection luminosities are given in the comoving frame. 

c Integrated 0.01 - 1 TeV y-ray luminosity of the model SED. Basically, this coincides with the observed value. 
d Observed total neutrino luminosity, i.e. (v^ + v^) + (v e + v e ). 
e Yyy is defined in equation G3. 

f This is the value derived without including in our modelling the upper limit in hard X-rays (see Sect. 14.13) . For this reason, it 
should be considered only as an upper limit. 

g Estimate of the pzr optical depth (or efficiency). We define it as f pn — 8 Z/ TeV /Z/ in j, where Z/ TeV is the integrated 0.01 - 1 TeV y- 
ray luminosity as measured in the comoving frame. The values are only an upper limit because of the assumption that the observed 
y-ray luminosity is totally explained by p/r interactions. 


the figures that follow, we will display, for clarity reasons, the n° 
component only before its attenuation by the internal synchrotron 
and EBL photons. 

lAbdo et alJ l l201ll) . whose data we use here, have also pre¬ 
sented a hadronic fit to the SED. Thus, a qualitative comparison 
between the two models is worthwhile. Both models are similar in 
that they require a compact emitting region and hard injection en¬ 
ergy spectra of protons and primary electrons. Morevoer, the y-ray 
emission from MeV to TeV energies, in both models, has a sig¬ 
nificant contribution from the hadronic component. However, the 
models differ in several aspects because of: (i) the differences in 
the adopted values for the magnetic field strength and maximum 
proton energy; (ii) the Bethe-Heitler process, which acts as an in¬ 
jection mechanism of relativistic pairs. Because of the s trong mag¬ 
netic field and large y p>max used in lAbdo et all d201 ill , the y-ray 
emission from Mrk 421 is mainly explained as proton and muon 
synchrotron radiation. In our model, however, these contributions 
are not important, and the y-rays in the Fermi/LAT and MAGIC 
energy ranges are explained mainly by synchrotron emission from 
pzr pairs (see DPM14). Moreover, the emission from hard X-rays 
up to sub-MeV energies is attributed to different processes. In our 
model, it is the result of synchrotron e missio n from B ethe- Heitler 
pairs, whereas in the hadronic model of lAbdo et alJ d201 lh . this is 


explained as synchrotron emission from pion-induced cascades. It 
is important to note that the B ethe-H eitler process was not included 
in the hadronic model of lAbdo et al.1 d20 1 ill . 

As far as the neutrino emission is concerned, we find that the 
model-predicted neutrino flux is below but close to the 1 a error bar 
of the derived value for the associated neutrino 9 (PR 14). This prox¬ 
imity is noteworthy and suggests that the proposed leptohadronic 
model for Mrk 421 can be confirmed or disputed in the near future, 
as IceCube collects more data. Even in the case of a future rejection 
of the proposed model for Mrk 421, there will still be room for lep¬ 
tohadronic models that operate in a different regime of the param¬ 
eter space. For example, in the LHs model, the neutrino spectrum 
is expected to peak at higher energies and to be less luminous than 
the one derived here, because of the higher values of y p , max and B 
(see also in DPMI4); this would also hold for the hadronic scenario 
used in lAbdo et ail fcOllh . although no discussion of neutrinos is 
presented therein. For such leptohadronic variants, ~PeV neutrino 
observations cannot be model constraining. 

4.1.2 PG1553+113 

The hybrid SED of blazar PG 1553+113 is shown in Fig. [4] The 
contribution of different emission processes to the total SED is de- 


©2013 RAS, MNRAS 000, [Tp?? 

































8 


leptonic + pe - 

leptonic + pn - 


pe — 

Pit - 

leptonic 

all processes - 


to — 

Ve+v„ . 

1C neutrino 1^ i— 1 —' 
KVA (min/max) • 
Swift/UVOT • 

logs (eV) 


Swift/XRT 
Swift/XRT-high 
RXTE 2008 
1FGL (2008-2009) 
2FGLJC 2008-2010) 
MAGIC (2005-2009) 


-4 -2 0 2 4 6 8 10 12 14 16 



Figure 4 . SEP of blaza r PG 1553+113 for the period 2005-2009 and the neutrino flux for the corresponding IceCube event (ID 17). The data are obtained 
from lAleksic et al]l2012l . Grey diamonds correspond to the KVA minimum and maximum fluxes and orange circles are the Swift/UVOT observations. Light 
and dark red diamonds are observations from Swift/XRT in 2005 corresponding to an intermediate and high flux level, respectively. With light green circles 
(from top to bottom) are plotted the average RXTE flux and the average 15-150 keV Swift/BAT flux. Green circles represent the average Fermi/LAT spectrum 
(August 2008- February 2009) and black circles the MAGIC observations over the periods 2005-2006 and 2007-2009. 


picted with different types of lines, which are explained in the leg¬ 
end above the figure. The same applies to all figures that follow. We 
did not do the same for the case of Mrk 421, since we quoted re¬ 
sults originally presented elsewhere, where a detailed description 
of the various processes contributing to the overall SED can be 
found. In any case, it was the prediction of the model about the neu¬ 
trino flux from Mrk 421 we wanted to emphasize and not the spec¬ 
tral modelling of the SED. A few things about the hybrid SED are 
worth mentioning. It is the SSC emission of primary electrons (blue 
dotted line) that mainly contributes to the GeV-TeV y-ray regime, 
whereas the p;r component is hidden below it (magenta dash-dotted 
line). Notice also that the other component of photohadronic origin, 
namely the synchrotron emission from Bethe-Heitler pairs (ma¬ 
genta dashed line), has approximately the same peak flux as the p n 
component. The fact that the SSC emission is favoured with respect 
to the emission from photohadronic processes has to do with the 
particular choice of parameter values (Table [2}. Besides the weak 
magnetic field (6 = 0.05 G), the combination of a large source 
(R = 2 x 10 17 cm) and a relatively high Doppler factor (t> = 30), de¬ 
creases the compactness of the primary synchrotron photon field, 
which is the target for the photohadronic interactions, and even¬ 
tually leads to a suppression of the proton cooling due to photo¬ 
hadronic processes and of the respective emission signatures (for 
more details, see Sect. 3.3 in PM15). 

In our framework, the total neutrino luminosity is expected to 
be of the same order of magnitude as the synchrotron luminosity 
of p n pairs that is emitted as high-energy y-rays (see also PM 15). 
This is clearly shown in Fig. [4] where the peak luminosity of the 
p n component is similar to the one of the neutrino component (red 
dotted line). Taking into account that the p7r component has a small 
contribution to the y-ray emission, we expect also the peak neu¬ 
trino flux to be lower than the peak y-ray flux, which is verified 
numerically (see Fig. 0. 

The previous results may be quantified through the ratio Y vy , 


which is defined as 

y n = 7 —, ( 13 ) 

Ly.TcV 

where L v and L yTe v are the total neutrino and 0.01 - 1 TeV y-ray 
luminosities, respectively, as measured in the observer’s frame. The 
ratio Y vy expresses the link between the observed y-ray luminosity 
and the model-predicted neutrino luminosity from a blazar, and as 
such is as an important parameter of this study; this will become 
clear later in Sect. 14.31 

For PG 1553+113, in particular, we find Y^ ~ 0.1, which 
suggests that the p7r component is not the sole contributor to the 
observed y-ray flux. Since the low value of Y vy is a matter of pa¬ 
rameter values, it might make the reader wonder why we adopted 
the specific parameter values in the first place, and why we did 
not search for parameter sets that would favour photopion produc¬ 
tion. There are two profound reasons for our particular choice: the 
shape of the Fermi/LAT spectrum, and the hard X-ray observations 
with Swift/BAT. Although a higher proton injection luminosity or a 
more compact emitting region would lead to higher neutrino fluxes 
and could, thus, account for the neutrino event ID 17, it would 
inevitably raise the Bethe-Heitler component. Since the emission 
from pe pairs falls, in general, in the hard X-ray and soft gamma- 
ray energy bands, the Bethe-Heitler component would not only ex¬ 
ceed the BAT observations but also destroy the fit to the low-energy 
Fermi/LAT data. One must also take into account that the pn syn¬ 
chrotron spectrum is broader around its peak than the SSC one, 
which would make it difficult to explain the observed Fermi/LAT 
spectrum. 

Inspection of Table|2]reveals that PG 1553+113 is the second 
most energetically demanding case, with L' m . being approximately 
one order of magnitude higher than the other values. In terms of 
parameter values, such as the size R and the magnetic field, it is an 
“outlier” of the sample. This should not come as a surprise, since 
PG 1553+113 is the most distant BL Lac of our sample with an 
estimated redshift z ~ 0.4 (see Sect.O, and its y-ray flux is approx- 
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imately equal to that of Mrk421. PG 1553+113 is, therefore, an in¬ 
trinsically more luminous source than Mrk 421 and requires an ac¬ 
cordingly higher proton injection luminosity as measured in the co¬ 
moving frame, unless a larger value for Doppler factor (6 > 30) was 
adopted. For a given R this would lower the required proton com¬ 
pactness, and thus Z/ mj . We attempted, however, to avoid higher 
Doppler factors than those listed in Table[2] since those would shift 
the neutrino peak to hundreds of PeV in energy (see equation 0). 

Even for the adopted parameters, the peak energy of the 
model-derived neutrino spectrum is ~ 10 PeV. If we combine this 
with the hardness of the spectrum we obtain a much lower flux than 
the observed one, at the energy of the associated neutrino event; 
still, this lies within the 3 cr error bars (see also Sect. 14.21 . 

In Fig. [4] we have also included the X-ray high flux state of 
2005 for illustrating reasons. An association of a higher X-ray flux 
with a higher neutrino flux would be intriguing, since it would bring 
the model-derived flux shown in Fig.[4]even closer to the observed 
one. However, such an association is not trivial in the absence of 
simultaneous observations in high-energies; neither GeV nor TeV 
data were ava ilable during the 2005 X-ray high state (see Fig. 1 
in lAleksic et al.ll2012l) . Consider for example a scenario where the 
increase in the X-ray flux is caused by an increase of the electron in¬ 
jection compactness. This would also result in a higher y-ray flux, 
since the emission from both the SSC and ptr components would 
increase. Finally, it would lead to a higher neutrino flux because of 
the higher synchrotron photon compactness. However, if the X-ray 
high flux state was not accompanied by an increase in the y-ray 
emission, the aforementioned scenario would not be viable. More 
than one model parameters should be changed, and, depending on 
their combination, the derived neutrino flux would also differ. Be¬ 
cause of the wide range of possibilities we do not, therefore, at¬ 
tempt to fit the high-state observed in the X-rays, nor to calculate 
the respective neutrino flux. 

4.1.3 1ES1011+496 

The hybrid SED is shown in Fig. [5] and the resemblance to the 
model SED of PG 1553+113 is obvious. The y-ray emission is 
mainly attributed to the SSC emission of primary electrons, since 
the flux of the p n component is only a small fraction of the ob¬ 
served y-ray one. The fact that the model neutrino peak flux is 
lower than the peak y-ray flux should come as no surprise after the 
discussion in Sect. l4.1.2l The model in Fig.[5]is mainly constrained 
by the simultaneous optical, soft X-ray and TeV data of the 2008 
MW campaign, since there are no significant spectral variations in 
the Fermi/LAT observations covering the period 2008-2010. As al¬ 
ready noted in PM 15, any information on the hard X-ray emission, 
e.g. e > 20 keV, would be even more constraining for the model. 
This becomes also evident from the cases of PG 1553+113 and 
H 2356-309 shown in Figs.[4]and[6] respectively. 

We searched, therefore, a posteriori for possible detections or 
upper limits in the hard X-ray band, even if these were not, strictly 
speaking, simultaneous with the 2008 data. For this, we used 
the fourth INTEGRAL/IBIS catalog in 20-40 keV lUbertini et alj 
120091) , that covers the period 2003-2008. The source has never been 
detected within this period, which resulted in the upper limit of 
~ 6.1 x 10- 12 erg cnr 2 s -1 shown with a grey arrow in Fig. [5] The 
model SED exceeds by a factor of ~ 3 the upper limit. This suggests 
that the photohadronic emission (both pe and p;r) should decrease, 
if we were to include the upper limit into our fitting procedure. 
Since the most straightforward way to decrease the photohadronic 
emission in order to accommodate the hard X-ray constraint is to 


assume a lower , the derived value for L p mj listed in Table [2] 
should be considered only as upper limit. A lower £ p j n j (keeping 
all other parameters unaltered) would also result in an accordingly 
lower neutrino flux and, thus, lower Yyy. Given that the model neu¬ 
trino flux at the same energy bin with neutrino event 9 is already 
below the observed one by approximately two orders of magnitude 
(see Figs. l5land lT0t . we argue that the hard X-ray constraint makes 
the connection of 1ES 1011+496 with the neutrino event 9 even 
less plausible. 


4.1.4 H 2356-309 

H 2356-309 is one of the sources in the sample that has also been 
detected in TeV energies by H.E.S.S. and is the most probable 
counterpart for the IceCube neutrino with ID 10. Our model SED 
and neutrino flux are shown in Fig. [0 In contrast to PG 1553+113, 
the Fermi/LAT spectrum can be approximated as vF v oc y° within 
the error bars and the y-ray flux at ~ 0.1 TeV is — 1 order of magni¬ 
tude below the peak flux of the low-energy hump. The properties of 
the observed SED of H 2356-309 allows us, therefore, to search for 
parameters that favour the ptr emission and suppress the contribu¬ 
tion of the SSC component, contrary to the cases of PG 1553+113 
and 1ES 1011+496. 

By adopting R = 3 x 10 15 cm, which is smaller by two and one 
orders of magnitude than the source’s radius used for PG 1553+113 
and 1ES 1011+496, respectively, we increase the photon compact¬ 
ness of the source. Moreover, the choice of a stronger magnetic 
field (B = 5 G) results in higher magnetic compactness (Zb) than 
before, where /g oc B 2 /R and is a measure of the synchrotron cool¬ 
ing rate. Thus, the synchrotron emission by secondary pairs from pe 
and p7r processes is enhanced compared to the previous two cases. 
Figure [6] displays the above. It is, indeed, the emission from the 
p;r pairs that mainly contributes to the observed y-ray flux, while 
the SSC component is ~ 1.5 orders of magnitude less luminous. 
Moreover, the “Bethe-Heitler bump” becomes a prominent feature 
of the SED. Both results suggest a higher efficiency of both photo¬ 
hadronic processes compared to the previous two cases. 

The neutrino spectral shape is similar to the previous cases and 
peaks at ~ 10 PeV. However, the ratio Y v . Y increases to 0.8, an eight¬ 
fold increase from the case of PG 1553+113, which reflects the dif¬ 
ferent efficiencies of photopion production in these cases. Although 
Yyy is close to unity, the peak neutrino peak flux is higher than the 
peak y-ray flux because of the different spectral shapes of the two 
components. At the energy of the detected neutrino (~ 0.1 PeV) the 
model-derived flux is low, but lies within the 3 cr error bars. This 
discrepancy is similar to the one we found in previous cases, such 
as PG 1553+113, although the photopion production efficiency and 
Y^ are both higher in the case of H 2356-309. At first sight, this 
result seems contradictory. However, one has to take also into con¬ 
sideration the differences in the ratios of the observed y-ray and 
neutrino fluxes. For example, in blazars PG 1553+113 and H 2356- 
309, this ratio is close to and much less than unity, respectively. 
Thus, in the case of H 2356-309, the higher Yyy is compensated by 
the lower ratio of observed y-ray and neutrino luminosities. 


4.1.5 1H1914-194 

Blazar 1H 1914-194 is one of the two sources that are potential TeV 
emitters and has the poorest MW coverage among the BL Lacs of 
the sample (for details, see Sect. [3j. Because of the inhomogene¬ 
ity (in time domain) of the available MW observations, we present 
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Figure 5. SED of blazar 1ES 1011+496 for the period 2008-2 010 and the neutrino flux for the corresponding IceCube event (ID 9). Optical, soft X-ray and 
TeV data were obtained during the first MW campaign in 2008 jReinthal et alJ2012h . whereas GeV observations are from the first (1FGL) and second (2FGL) 
Fermi/LAT catalogs (quasi-simultaneous). The first TeV detection is over-plotted with grey points for comparison reasons, as it corresponds to a high flux state; 
the respect ive data points were obtained from ASDC. The grey arrow denotes the upper limit to the 20 keV-40 keV X-rays from the fourth INTEGRAL/IBIS 
catologue iUbertini et all20091) . 
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Figure 6. SED of blazar H 2356-309 for the period 2004-2010 and the neutrino flux for the corresponding IceCube event (ID 10). The data are obtained from 
ASI Science Data Center using the SED builder tool. All types of symbols and lines are explained in the legend above the plot. 


only an indicative fit, where we attempted to describe at best the 
Fermi/LAT data. Being the least constrained case we chose param¬ 
eters that ensure high efficiency in photopion production, i.e. we 
present the most optimistic, for neutrino emission, scenario. We 
caution the reader that contemporaneous data in optical, X-rays and 
GeV y-rays may strongly constrain the model (see e.g. Mrk 421 
and PG 1553+113) and lead to different parameter values from 
those listed in Table [2] The neutrino flux shown in Fig. [7] should 
be therefore considered as an upper limit. We also note that the 
non-detection of TeV y-rays is not as important for the predicted 
neutrino flux as the simultaneous coverage in the aforementioned 
energy bands. For example, the peak neutrino energy is solely de¬ 


termined by the peak frequency of the low energy component and 
the value of the Doppler factor (see equation 0), while the neu¬ 
trino flux is related to the sub-TeV y-ray flux (see Figs. am and 
Fig. 8 in PM15). 

The SED we obtained is the most complex one in terms of 
the number of emission components that appear in the MW spec¬ 
trum. Apart from the very low frequencies (below UV) where the 
primary leptonic component contributes, the spectrum is comprised 
roughly of: proton synchrotron radiation (between UV and ~keV), 
synchrotron emission from Bethe-Heitler pairs (between ~keV and 
~100 MeV) and synchrotron emission (< 1 TeV) from secondary 
pairs produced through pion decay and yy absorption of VHE n° 
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Figure 7. SED of blazar H 1914-194 and the neutrino flux for the corresponding IceCube event (ID 22). The data are obtained from ASI Science Data Center 
using the SED builder tool. All types of symbols and lines are explained in the legend above the plot. 
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Figure 8. SED of blazar 1 RXS J054357.3-553206 and neutrino flux for the corresponding IceCube event (ID 19). The data are obtained from ASI Science 
Data Center using the SED builder tool. All types of symbols and lines are explained in the legend above the plot. 


y-rays (see equation (14) in lPetropoulou & Mastichiadisll2012h . In 
fact the role of yy absorption is non-negligible in this case, since 
the unattenuated 7r° component has a flux higher by a factor of ~ 2 
than the neutrino one (see dash-dotted gr ey curve). The attenua - 
tion of VHE y-rays initiates an EM cascade ((Mannheim et all 19911) 
that transfers energy from the PeV regime to lower energies. Our 
numerical approach does not allow us to isolate multiple photon 
generations produced in the electromagnetic cascade. Thus, strictly 
speaking, the pe and pn components shown in Fig. [7] are the re¬ 
sult of synchrotron radiation not only from pairs produced through 
photohadronic interactions but also from pairs produced in the EM 
cascade. 

In this case, we find Yyy - 2, which is the highest value 
amongst the sample (see Tabled and close to the maximum value 
(T“ ax = 3) that is allowed in this theoretical framework (see Sect. 


3.2 in PM15). The model-predicted neutrino flux is close to the ob¬ 
served flux (within the 1 cr error bars), thus strengthening the asso¬ 
ciation of neutrino event 22 with this source. We have to note, how¬ 
ever, that the model-predicted neutrino flux might decrease, should 
more constraints on the SED be available in the future. 

4.1.6 1 RXS J054357.3-553206 

1 RXS J054357.3-553206 is the second HBL source that has not 
yet been detected in TeV energies and the most probable counter¬ 
part of IceCube event 19 according to PR14. The hybrid SED ob¬ 
tained with our model is shown in Fig. H] and shares many features 
with the model SED of PG 1553+113. According to the discussion 
in Sect. 14.1.21 we can infer that the observed y-ray emission origi¬ 
nates only partially from photohadronic processes, since the model 
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derived neutrino peak flux is lower than the y-ray one. In fact, we 
find that Y n = 0.1, same as in the case of PG 1553+113. However, 
because the flux corresponding to neutrino event ID 19 is higher 
than the y-ray one, the model-predicted flux (at the same energy 
bin as event 19) is more than 3 cr off from the observed value, and 
the model prediction falls short of explaining the observed flux (see 
also Fig.ITTfli 

In contrast to PG 1553+113 the SED of blazar 
1 RXS 1054357.3-553206 is less constrained both in (sub)MeV 
and TeV energies. This, along with the large discrepancy between 
the model and observed neutrino fluxes, makes 1 RXS J054357.3- 
553206 the best case for studying the effects of a different 
parameter set on the neutrino flux. We search, therefore, for 
parameters that increase the contribution of the p7r component 
to the observed y-ray emission and at the same time lead to a 
higher neutrino flux. We do not focus, however, on the goodness 
of the SED fit but we rather try to roughly describe the observed 
spectrum. The aim of what follows is (i) to demonstrate the flex¬ 
ibility of the leptohadronic model in describing the MW photon 
spectra and (ii) to find parameters that maximize the predicted 
neutrino flux from the particular blazar, albeit at the cost of a worse 
description of the SED. 

Our results are summarized in Fig. [9] where we plot the hy¬ 
brid SED obtained for a different parameter set presented in Ta¬ 
ble]^] From this point on, we will refer to it as Model 2. For a better 
comparison, we over-plotted with dashed lines the photon and neu¬ 
trino spectra that were obtained for the parameters listed in Table [2] 
(Model 1). In Model 2 the y-ray emission is mainly attributed to 
the pn component and the neutrino flux is comparable to the ob¬ 
served y-ray one; in particular, we find Yyy ^ 3. In our framework, 
this is the most optimistic scenario and the derived neutrino flux 
can be considered as an upper limit on the flux expected from this 
source (see also discussion in Sect. 14. 1 .5l >. However, the description 
of the SED is not as good as in our baseline model. In particular, 
the model-derived spectrum (i) cannot explain the highest energy 
Fermi/LAT data, (ii) is less hard (in vF v units) at GeV energies and 
more broad than the observed one - see also Sect. 14.1.21 and (iii) 
is marginally above the highest flux observation in X-rays. This 
excess is caused by the synchrotron emission from Bethe-Heitler 
pairs, which appears in the hard X-ray/sub-MeV energy range, and 
it is even more luminous than the pn component. It is also note¬ 
worthy that the Bethe-Heitler emission in Model 2 plays a crucial 
role in the neutrino production: if we were to neglect pair injection 
from the Bethe-Heitler process, the neutrino flux in Model 2 would 
be lower by a factor of ~ 10. The detection of a luminous compo¬ 
nent in the 20 keV-20 MeV energy range could, therefore, serve as 
an indirect probe of high-energy neutrino emission from a blazar 
source within the leptohadronic scenario, as recently proposed in 
PM15. 

4.2 Neutrino emission 

The model-derived electron and muon v + v spectra for the six 
blazars of the sample are summarized in Fig. [TO] For a better com¬ 
parison, we over-plotted the fluxes that correspond to the respective 
neutrino events along with the Poissonian lcr error bars calculated 
for a single event (PR14). We note that re spective 3 cr error bars 
correspond to -2.9, +0.9 dex dGehrelsfl986l) . 

We focus first on the neutrino event 9, which according to the 
model-independent analysis of PR 14, has two plausible astrophys- 
ical counterparts; namely blazars Mrk 421 and 1ES 1011+496. For 
Mrk 421, in particular, we show the neutrino spectrum depicted in 


Model 1: photons . 1C neutrino 19 i—•—i 1FGL (2008-2009) 

Model 2: photons - Swift XRT i— ♦—i 2FGL (2008-2010 

Model 1: v e +v - WISE 2010 + radio i—•—i 2FGL_lc (2008-2010) 

Mortal O- \t t\> vhAhA /onnc nr\r\~r nr\-tr\\ 
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logs (eV) 
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Figure 9. SED of blazar 1 RXS J054357.3-553206 and neutrino flux for 
the corresponding IceCube event (ID 19). The hybrid model SEDs obtained 
for two different parameter sets presented in Tables[2]and[3]are respectively 
shown with dashed and solid lines. 


the right panel of Fig. [3] According to the discussion in Sect. [5] 
the neutrino spectrum for 1ES 1011+496 (dashed line) should be 
considered as an upper limit, since inclusion of additional con¬ 
straints in the SED fitting would lower the neutrino flux by at least 
a factor of ~ 3. Thus, our results strongly favour Mrk 421 against 
1ES 1011+496. As already discussed in the respective paragraphs 
of Sect. [4] the differences between the neutrino fluxes originate 
from the differences in their SEDs. Thus, the case of neutrino ID 9 
reveals in the best way how detailed information from the photon 
emission may be used to lift possible degeneracies between multi¬ 
ple astrophysical counterparts. 

1 RXS 1054357.3-553206 was selected by PR14 as the as¬ 
trophysical counterpart of neutrino event 19. Figure [TO] shows 
the respective neutrino spectrum (yellow line) we obtained for 
the best-fit model of its SED. The model-derived flux at the en¬ 
ergy of ~ 0.2 PeV lies below the 3 cr error bars. In this regard, 
1 RXS 1054357.3-553206 may be excluded at the present time 
from being the astrophysical counterpart of its associated neutrino 
event, namely event 19. We note, however, that a higher neutrino 
flux (within the 3 cr error bars) can be obtained at the cost of a 
worse description of the SED (see Sect. [81. 

In all other cases, the model-deriveqj neutrino flux at the en¬ 
ergy bin of the detected neutrino is below the 1 cr error bars, but still 
within the 3 cr error bars. Thus, strictly speaking, the association of 
these sources cannot be excluded at the present time. In particular, 
Mrk 421 and 1H 1914-194 are the two cases with the smallest dis¬ 
crepancy between the observed and model-derived fluxes and, in 
this regard, the most interesting, because their association with the 
respective IceCube events can be either verified or disputed in the 
near future. 

Besides the comparison between the model and the observed 
neutrino fluxes. Fig. [TO] also demonstrates the similarity of the 
model neutrino spectra, despite the fact that they were obtained 
for BL Lacs with different properties. Our results suggest that the 
shape of the neutrino spectrum depends only weakly on the y-ray 
luminosity of the source, which varies approximately two orders of 
magnitude among the members of the sample (see Table[2j- In par- 


7 We refer to the neutrino emission calculated for best-fit SED models. 
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1ES 1011+496 (ID 9) 
PG 1553+113 (ID 17) 


H 2356-309 (ID 10) 
1H 1914-194 (ID 22) 
1RXS J05435-5532 (ID 19) 



Figure 10. Comparison of the model (lines) and the observed (circles) neutrino fluxes as defined in PR 14 for the six BL Lacs of the sample. The Poissonian 
1 cr e rror bars for each event, which are adopted by PR 14, are also shown. We note that the respective 3 cr error bars correspond to -2.9, +0.9 dex iGehrelsI 
Il98+ 


ticular, the neutrino spectrum may be described as a power-law, i.e. 
E;dN v /dE v cc E 1 ’ with a ~ 1.1 - 1.3. 

Before closing this section, we must point out that the present 
theoretical framework allows us to establish a connection between 
the observed y-ray luminosity and the predicted total neutrino lu¬ 
minosity. In fact, we showed that these are simply proportional as 
L v = YyyLyj c v, with Yyy in the range 0.1-2 (Table 13. Although 
there is a maximum value predicted by the model, i.e. Yyy < 3 (see 
also PM 15), there is no restriction on the lowest possible value of 
this ratio. Values of Yyy <sc 1 will only imply that the leptohadronic 
model simplifies into a pure leptonic one, with the SSC emis¬ 
sion of primary electrons being responsible for the observed emis¬ 
sion. This is illustrated best through the cases of PG 1553+113, 
1ES 1011+496, and 1 RXS J054357.3-553206. 

4.3 The ratio Yyy and the energetics 

In this section we examine in more detail the results presented in 
Table[2] We attempt to extract information about the dependence of 
the ratio Y vy , which is a crucial derivable quantity of our model, on 
observables, such as the photon index in the Fermi energy range and 
the y-ray luminosity L y _ TeV . Moreover, using the parameter values 
of Table[2] we calculate the jet power for each source and comment 
on the energetics. The results that follow, particularly those related 
to trends and correlations, should be considered with caution be¬ 
cause of the limited size of the sample. Although these cannot be 
directly applied to a wider sample of BL Lacs, they are useful in 
that they provide a better understanding of the modelling results. 

In the top panel of Fig. we plot th e photon index in the 
Fermi/LAT energy range iAbdo et aPl201(l) versus the ratio Yyy. 
Sources with ratio values above and below the mean value of the 
sample Y vy = 0.8 are plotted with black and red symbols, re¬ 
spectively. We remind the reader that for 1ES 1011+496 the ratio 
Yyy = 0.5 was derived without including in our analysis the up¬ 
per limit to the hard X-rays. As discussed in Sect. [5] inclusion of 
the hard X-ray constraint would lower the current ratio, at least by 


Table 3. Parameter values used for the fit presented in Table [2] (Model 1) 
and in the fiducial model described in text, which maximizes the neutrino 
tlux (Model 2) from blazar 1 RXS J054357.3-553206. The symbols have 
the same meaning as in Table[2] 


Parameter 

Model 1 

Model 2 

B 

0.1 

1 

R (cm) 

3 X 10 16 

3 x 10 15 

<5 

31 

50 

ye,min 

1 

1 

Te,max 

1.2 x 10 5 

2.5 x 10 4 

s e 

1.7 

1.7 

^e,inj 

2.5 x 10- 5 

1.6 x 10~ 5 

yp,min 

1 

1 

Tp,max 

1.2 x 10 7 

1.2 x 10 6 

i P 

2.0 

1.9 

f ■ ■ 
l/ P>inj 

8 x 10~ 4 

1.6 x 10~ 2 

Ty.TcV (erg/s) 

6.3 x 10 45 

1.6 x 10 45 

Ly (erg/s) 

6.3 x 10 44 

5 x 10 45 

Yyy 

0.1 

3.1 


a factor of ~ 3; this is illustrated with an arrow in both panels of 
Fig-EH There is no evident trend between the ratio Yyy and the pho¬ 
ton index in the GeV energy range, although in Sect. l4T.3 we used 
the hardness of the y-ray spectrum as an argument for choosing 
parameters that favoured SSC against pn radiation. These findings 
suggest that the Fermi photon index alone is not a good indicator 
of the model-derived ratio Y vy . However, we find a negative trend 
between the observed y-ray luminosity in the 0.01-1 TeV energy 
range (listed in Table[2} and the ratio Y vy (bottom panel in Fig. II 111 

8 One could argue that we were bound to find an anti-correlation, since we 
are plotting L y ,xeV against the ratio Y v} , (note that Z. r ,T e v = T,7 y L v ). This 
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Figure 11. Top panel: Plot of the photon index in the Fermi/LAT energy 
range versus the ratio Y vy . Sources with ratio values above and below the 
mean value Y yv = 0.8 are plotted with red and black symbols, respectively. 
The arrow stresses the fact that the ratio derived for blazar 1ES 1011+496 
is only an upper limit (see Sec t. [31 . The values of the photon index are 
adopted from the 1FGL catalog I.Abdo et al]|2010l) . Bottom panel: The ob¬ 
served (0.01-1) TeV y-ray luminosity plotted against the ratio Y vy (logarith¬ 
mic scale). The arrow and black (red) symbols have the same meaning as in 
the top panel. The dashed line is a linear fit (log L y> t c v = A 1 log Y vy + Ao) 
to the points with Aj = -1.1 ± 0.4, Aq = 45.2 ± 0.2. 


The size of our sample is limited, yet our results suggest that the 
contribution of the pn component to the blazar’s y-ray emission is 
smaller in sources that are more y-ray luminous. The dependence 
of Yyy on the y-ray luminosity is particularly useful if we attempt to 
calculate the BL Lac diffuse neutrino flux by extrapolating the find¬ 
ings of this sample to the wider class of BL Lacs; this will be part of 
a future study. Finally, we point out that the negative trend shown 
in Fig.QT]does not necessarily mean that less powerful BL Lacs in 
y-rays inject more power in relativistic protons, as it will become 
evident at the end of this section. 

The bottom panel of Fig.|TT]provides the first hints of a separa¬ 
tion between the sources of the sample with ratios above and below 
the average value. Flowever, this becomes clear in Fig. EQ where 
we plot the proton (filled symbols) and electron (open symbols) 


would be, indeed, correct, if there was an a priori knowledge of the model- 
derived neutrino luminosity. In fact, there was no guarantee that the neutrino 
luminosity obtained for all six BL Lacs would be approximately constant 
(see Table [2), which is an interesting outcome of the present analysis. 


Figure 12. Logarithmic plot of the proton (filled symbols) and electron 
(open symbols) injection compactnesses as a function of the magnetic com¬ 
pactness (for the definition, see text) for the six BL Lacs of the sample. 
Red (black) colored symbols correspond to ratio values Yyy above (below) 
the average value derived from the sample, i.e. Yyy = 0.8. The ^ Pi i n j value for 
1ES 1011+496 is an upper limit, shown with an arrow. The dashed arrow 
demonstrates the transition from SSC to pn dominated y-ray emission. 



lo g L p,inj( er g /sec ) 


Figure 13. L' versus L' in logarithmic scale. The color coding is the 

° e.inj_ p,inj ° ° 

same as in Fig. 1121 The solid and dashed lines have slopes equal to one and 
two, respectively, and are plotted for guiding the eye. 


injection compactnesses against the magnetic compactness defined 
as £b = crjRB 2 /8mn t c 2 . Red and black colored symbols are used 
to denote values of Y vy above and below the mean value, respec¬ 
tively. Figurefl2l reveals the presence of two sub-groups within the 
sample. As increases, the synchrotron cooling of both primary 
and secondary electrons is enhanced. On the one hand, this leads to 
an increase of the synchrotron emission from pairs produced by p;r 
interactions that emit in the y-ray regime, while on the other hand, 
it suppresses the SSC emission of primary electrons. Thus, as 
increases, the contribution of the pn component to the y-ray emis¬ 
sion increases against the SSC one. This is finally reflected in an 
increase of the ratio Y vy . 

In Fig.|T3]we plot L e mj as a function of L p mj . The color coding 
is the same as in Fig. ED and the arrow stresses the fact that the 
derived value for 1ES 1011+496 is an upper limit (see Sect. 14.1.3k 
The injection luminosities of primary particles, as measured in the 
comoving frame, follow a more than linear relation, i.e. L e mj oc 
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(L p m ^ with 1 < x < 2. More important, though, is that sources 
with Yyy > 0.8 lie at the lower left corner of the plot, i.e. they are 
described by low proton and electron luminosities (in the comoving 
frame) with respect to the BL Lacs with smaller F vy . 

The jet power is a useful quantity when it comes to the ener¬ 
getics of the model. In the one-zone leptohadronic model, the en¬ 
ergy densities of both relativistic electrons (n') and protons ( u ' p ), as 
well as the energy density of the magnetic field (u B ) and radiation 
(nj), contribute to the jet power, which may be written as 

Lj =: 7tR 2 8 2 pc (u' e + u p + u' B + n'), (14) 

where the energy densities of the various components are mea¬ 
sured at steady state and in the comoving frame, 6 ~ T, and 
y3 = yT— 1/T 2 g 1. For the assumptions u sed to derive equa¬ 
tion ( 1141 . see e.g. ICelotti & GhisellinJ d2008l) . The jet luminosity 
as well as u[ are listed, for reference, in Table 0] For all sources we 
find high jet luminosities, and emission regions far from equipar- 
tition u' p /u' B » 1, with the proton component being the domi- 
nant one. The s e resu l ts are typical of lep tohadronic models (e.g. 
iBdttcher et alJ d2013h : ICerruti etaD d20 1 4l ) ) and often constitute a 
point for criticism. However, our analysis highlighted that sources 
with Yyy ~ 1, namely with significant contribution of the photo- 
hadronic component to the y-ray emission, are described by less 
extreme, i.e. < 10 49 erg/s, jet powers and have lower u p /u' B ratios; 
this was not expected beforehand. For the rest of the sources, i.e. 
those with T vy <K 1, the derived jet luminosities exceed 10 49 erg/s, 
mainly because of the larger R compared to the other sources (see 
Table 0 ). 


5 DISCUSSION 

In the present paper we have applied a one-zone leptohadronic 
model to six BL Lacs, which have been recently selected as the 
probable astroph ysical counterparts (PR 14) of five of the IceCube 
neutrino events dAartsen et alJl2014D . and calculated the neutrino 
flux in each case. In contrast to other studies, where the neutrino 
emission is calculated under the assumption of a generic blazar 
SED, here we obtained the neutrino signal from BL Lacs by fit¬ 
ting their individual SEDs. Thus, for each source, the photon and 
particle distributions used to derive the neutrino emission have to 
be determined individually, so as to give the observed SEDs. 

For the SED modelling, we have adopted one of the possible 
leptohadronic model variants that allows us to directly associate the 
y-ray blazar emission with a high-energy (2-20 Pe V) neutrino sig¬ 
nal. In particular, the low-energy hump of the SED is explained by 
synchrotron radiation of primary relativistic electrons, whereas the 
observed high-energy (GeV-TeV) emission is the combined result 
of synchrotron self-Compton emission (from primary electrons) 
and of synchrotron radiation from secondary pairs. These are pro¬ 
duced in a variety of ways, such as in Bethe-Heitler and yy pair 
production, and through the decays of charged pions, which them¬ 
selves are the by-product of p n interactions of co-accelerated pro¬ 
tons with the internally produced synchrotron photons. Secondary 
pairs produced in different processes have, in general, different en¬ 
ergy distributions. Thus, in our framework, it is mainly the syn¬ 
chrotron emission of pn pairs that falls in the y-ray regime, and in 
combination with the SSC radiation is used to explain the observed 
blazar emission in y-rays. 

BL Lacs are known for their variability across the electromag¬ 
netic spectrum and, because of this, we used SEDs that are com¬ 
prised of nearly simultaneous MW observations. For those sources 


that were targets of recent MW observing campaigns we have 
adopted the data from the respective publications. In any other case, 
we have compiled the SEDs using the SED builder tool of the ASI 
Science Data Centre (ASDC). Only for one source in the sample 
(1H 1914-194), which had the poorest simultaneous MW cover¬ 
age, had we to construct an SED using non-simultaneous data. We 
used the numerical code described in DMPR12 in order to calcu¬ 
late steady-state photon spectra that were later applied to the MW 
observations. Because of the multi-parameter nature of the problem 
(11-13 free parameters) and of the long numerical simulation time 
required for every parameter set, we have not performed a blind 
search of the available parameter space. To obtain a first estimate 
of the required parameter set, we used instead the analytical expres¬ 
sions presented in PM15. We then performed a series of numerical 
simulations using parameter values lying close to the initial param¬ 
eter set, until a reasonably good fit to each SED was obtained. 

Although the model-derived neutrino emission from the six 
BL Lacs that were selected as counterparts of the IceCube events 
was the main focus of this paper, our study highlighted also the 
flexibility of the leptohadronic model and its success in describing 
the blazar SEDs. This is particularly important, if one takes into 
account the following: (i) all sources in the sample were fitted using 
reasonable parameter values; (ii) these BL Lacs were not a priori 
selected for modelling purposes; and (iii) besides the fact that all 
sources in the sample were HBLs, they differed in many aspects, 
e.g. in redshift, in y-ray luminosity and spectra. 

The only source that may be excluded, at the present time, 
from being the astrophysical counterpart of the detected neutrino 
(event 19) is 1 RXS J054357.3-553206. In all other cases, the 
model-predicted neutrino flux is below the 1 a but still within the 
3 cr error bars of the respective IceCube neutrino events. In partic¬ 
ular, for blazars Mrk 421 and 1H 1914-194 the model-derived neu¬ 
trino fluxes at the same energy with the respective neutrino events 
(ID 9 and 22, respectively) were found to be very close to the low 
1 cr error bars (Fig. [H. Although our results cannot be considered 
a proof of a firm association between these BL Lacs and the respec¬ 
tive neutrinos, they can be verified or disputed in the near future, as 
IceCube collects more data and the observed fluxes become more 
constraining for the model. 

Among the blazars of the sample, Mrk 421 is a special case. 
Its neutrino emission was calculated in a previous study (DPM14), 
at a time when there was no hint of a possible association between 
Mrk 421 and the IceCube events. DPM14 found that the neutrino 
flux from Mrk 421 was close to the current IceCube (IC-40) sensi¬ 
tivity limit for this source (left panel in Fig.[3]l. Of particular interest 
are the following: 

• the prediction of DPM14 was found a posteriori to be close to 
the observed flux related to neutrino ID 9; 

• the same parameter set, besides the Doppler factor, used in 
DPM14 to fit the optical/X-ray/TeV data of 2001, has been success¬ 
fully (right panel in Fig.[3]l appl ied t o the more complete dataset of 
the 2009 campaign lAbdo et ahkoioll . 

Apart from a direct comparison of the model and observed 
neutrino fluxes, we have shown that we can establish a direct con¬ 
nection between the observed y-ray emission and the putative neu¬ 
trino emission from a BL Lac source. We have quantified this re¬ 
lation by introducing the ratio of the total neutrino to the 0.01 - 1 
TeV luminosity (Tyy). In our model, for a given observed y-ray lu¬ 
minosity that is purely explained by the synchrotron emission of p n 
pairs, i.e. when the SSC contribution is negligible, there is an up¬ 
per limit to the neutrino luminosity. It can be shown that the limit 
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Table 4. Jet luminosity as defined in equation 1141 and (comoving) energy densities of the proton, electron, radiation and magnetic components. 


Source 

Lj (erg/s) 

wj, (erg/cm 3 ) 

«e (erg/cm 3 ) 

u[ (erg/cm 3 ) 

u' B (erg/cm 3 ) 

Mrk 421 

2 X 10 48 

4 x 10 3 

4 x 10~ 3 

4.9 x 10^ 2 

9.9 x 10 _1 

PG 1553+113 

3.7 x 10 49 

1.1 x 10 

4 x 10~ 4 

2.3 x 10~ 4 

9.9 x 10~ 5 

1ES 1011+496 

8.3 x 10 49 

9x 10 2 

6.6 x 10~ 3 

5.3 x 10- 4 

4 x 10~ 4 

H 2356-309 

4.4 x 10 48 

5.7 x 10 3 

1.4 x 10~ 3 

1.6 x 10^ 2 

9.9 x 10“' 

1H 1914-194 

8 x 10 48 

2.9 x 10 4 

7.4 x 10~ 2 

8.2 x 10^ 2 

9.9 x 10“' 

1 RXS J054357.3-553206 

1.5 x 10 49 

1.8 x 10 2 

3.3 x 10~ 3 

7 x 10~ 4 

4 x 10~ 4 


is ~ 3L y ,Tev (see PM15). If, however, the SSC component is domi¬ 
nant in the y-ray regime, we expect Yyy —» e <K 1; in this case, our 
model simplifies into a leptonic one. Thus, the spectrum of possi¬ 
ble values is [e <sc 1,3] and corresponds to different contributions 
of the p n and SSC components to the y-ray emission. 

For the six blazars we have found ratios in the range 0.1 < 
Yyy < 2, with a mean value Yyy — 0.8. We have identified two sub¬ 
groups in the sample: blazars Mrk 421, H 2356-309 and 1H 1914- 
194 with Yyy Yyy and the rest with ratios below the average 
value (Table |2j. Despite the small size of our sample, we have in¬ 
vestigated the dependence of this ratio on observables of blazar 
emission. No hints of correlation with the photon index in the 
Fermi/LAT energy range were found (Fig. |TT] top panel). However, 
we showed that there is a negative trend between L 7 ,tcv and Y vy , i.e. 
more y -ray luminous BL Lacs have a smaller photohadronic con¬ 
tribution to their y -ray spectra (bottom panel in Fig. ITU . Although 
this could be interpreted as a lower injection luminosity of relativis¬ 
tic protons, we showed that this is not the reason (see Fig.ll3land 
Table0. Instead, the separation of sources according to Yyy can be 
understood in terms of the magnetic compactness. In Fig. ED we 
demonstrated the transition from SSC to p7r dominated y -ray emis¬ 
sion as the magnetic compactness gradually increases, which also 
translates to an increase of Y vy . Application of the model to a larger 
sample of BL Lacs is, however, necessary for a better understand¬ 
ing of the aforementioned trends. 

Our analysis showed that the shape of the neutrino spec¬ 
trum does not strongly depend on the y -ray luminosity. In all 
cases, the neutrino spectra are described as dN v /dE v oc £" -2 with 
a — 1.1 — 1.3 (Fig. ED. The peak energy of the neutrino spec¬ 
trum, in the present context, depends only on one model parameter, 
namely the Doppler factor, and two observables, i.e. the peak fre¬ 
quency (v s ) of the low-energy hump of the SED and the redshift 
of the source z. The dependence on the latter is not strong. Since 
all the sources of the sample are HBLs and we have obtained fits 
with similar values of 6, we should not expect large differences 
in the peak neutrino energy. We have verified this numerically, as 
demonstrated in Fig. m In addition to the similarity in the spectral 
shape, we have connected the expected neutrino luminosity from a 
BL Lac with its y-ray luminosity in the 0.01-1 TeV energy range. 
Thus, one could go one step further and calculate the diffuse neu¬ 
trino flux from BL Lacs, in order to see how it compares with the 
current IceCube detections and upper limits in the sub-EeV regime. 
Such a calculation requires an extrapolation of our findings to the 
wider class of BL Lacs, including low-synchrotron peaked blazars 
(LBLs), and the use of detailed simulation s of the blazar population 
of the type done by IPadovani & Giommil J2015I ). This will be the 
focus of a future study. 

Similarly to other leptohadronic models, such as the proton 
synchrotron model (e.g. lBottcher et all2013l) . the derived jet power 
is high, ranging from 10 48 erg/s up to a few 10 49 erg/s, which im¬ 


plies a very low radiative efficiency. Moreover, the emitting re¬ 
gion is far from equipartition, with the proton energy density being 
the do minant one - see Tables [4] in this work and IBottcher et alj 
d2013l) . The demanding energetics and the large departure from 
equipartition constitute the main point for criticism of our model, 
and of leptohadronic models in general. However, there are some 
important differences compared to the proton synchrotron model, 
which make this leptohadronic variant (LH;r) more flexible. On 
the one hand, strong magnetic fields (» 10 G) are not necessary, 
while, on the other hand, proton acceleration to ultra-high energies 
(y p ~ 10 9 - 10 10 ) is not required. On the contrary, the upper cutoff 
of the proton distribution should be just a few times higher than 
~ 3.5 x 10 7 (1 + z) -1 £iv:} 6 (see equation 0). In terms of the avail¬ 
able parameter space of leptohadronic models, the LHtt variant oc¬ 
cupies a sub-region, which corresponds, roughly speaking, to: low- 
to-moderate magnetic field strengths, moderate maximum proton 
energies, sub-pc sizes of the emission region (~ 10 15 - 10 16 cm) 
and high jet powers. 

A detailed calculation of the neutrino emission from blazars 
following the so-called “blazar sequence” ( GhisellinyetalJ J_998j) 
has been presented in tw o recent studies by Murase et alj 1 20141) 
and lDermer et alj d2014h . Our analysis showed, however, that rel¬ 
atively high neutrino luminosities (of the order of ~ 10 45 erg/s) 
can be obtained from classical BL Lacs, with Mrk 421 being the 
most notable case. We note also that the p7r production efficiency 
is low in all six blazars (see Table[2]i, and is in agreement with the 
expec ted values from BL Lacs (see equation (20) in lMurase et alj 
120141 and equation (28) in PM 15). The question is then whether 
an d how those r esults are compatible with each other. Figure 9 
in Murase et alj 1201 4l) shows that the typical neutrino luminos¬ 
ity of individual BL Lacs with no contamination by external pho¬ 
ton fields (two lower blue curves) are ~ 10 42 erg/s, which is ap¬ 
proximately 3 orders of magnitude below our typical values. The 
reason behind this discrepancy is the luminosity that is being in¬ 
jected into relativistic protons. In our analysis, we systematically 
used 2 to 3 orders of magnitud e higher proton injection lumi¬ 
nosities than iMurase et ali ( 120141) . This is mainly due to the fact 
that we try to model the observed blaza r y-ray emissi on in terms 
of photohadronic p rocesses, whereas in IMurase et alj d2014l) and 
Der mer et ~aD J2014) it is not clear what the fate of the photons pro¬ 
duced via photohadronic processes is, and why they do not appear 
in the observed SED. In any case, we compensate for the low num¬ 
ber density of internal radiation by pushing the proton injection 
luminosity to higher values. 

Another way to go around the problem of low photon number 
density in BL Lacs, which is related to low ne utrino pro duction effi¬ 
ciency, has been also recently presented in lTavecchio et al.l d2014l ). 
Their approach is different than ours, though. In their model, 
BL Lacs can still be efficient neutrino emitters if the number density 
of photons, which are the targets for p n interactions, is enhanced. 
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In order to achieve this. lTavecchio et all 120141) assume a structured 
BL Lac jet with a fast spine a nd slow layer i Ghisellini et alj|2005l ; 
iTavecchio & Ghisellinil2008ll . where the target field for p n interac¬ 
tions is the radiation field from the layer that appears boosted in the 
rest frame of the spine. 

The recent correlation analysis of IceCube neutrino events and 
y-ray sources by PR14 led to a somewhat unexpected result: instead 
of the more powerful FSRQs, it was BL Lacs that were among 
the most plausible astrophysical counterparts of 9 IceCube neutri¬ 
nos. In general, FSRQs are expected to be efficient PeV neutrino 
emitters, because of their higher y-ray luminosity, but more im¬ 
portantly, because of the dense external photon fields (broad line 
region and/or dusty torus). The absence of FSRQs from the list of 
probable astrophysical counterparts may be related to the “ener¬ 
getic criterion” applied by PR14. This was totally empirical and 
observationally driven, since it assumed only a tight photon - neu¬ 
trino connection; any other more sophisticated choice of a diag¬ 
nostic would be model-dependent. Thus, sources where a simple 
extrapolation of the y-ray flux to the PeV energies lied much be¬ 
low the observed neutrino flux (taking into account the uncertainty 
in the flux measurement) were discarded. Relaxing this constraint 
would certainly increase the number of candidates, which might in¬ 
clude FSRQs, but these would then depend on the specific model 
assumed. 

In principle, the leptohadronic variant we presented here can 
also be applicable to FSRQs, with the appropriate choice of pa¬ 
rameters. Because of the higher y-ray luminosities and lower peak 
energies of FSRQs with respect to HBLs, a straightforward way to 
apply our model to FSRQs would be to increase the proton injection 
luminosity, while lowering the maximum proton energy. Taking 
into account the presence of additional photons originating from 
the broad line region and/or the torus, the model-predicted (PeV) 
neutrino luminosity would be much higher than the one derived for 
the HBLs in our sample and possibly close to, or at, the limit of flux 
detected by IceCube. The application of the model to the SEDs of 
FSRQs, however, may not be as simple as described above. Higher 
L' . ■ leads not only to increased emission from p7r interactions but 

p,mj J r 

also from Bethe-Heitler interactions. Given its spectral shape and 
the typical energy range that emerges (see Figs. I3I9L it would be 
difficult to reconcile its emission with the flat X-ray spectra of FS¬ 
RQs. However, more concrete answers can be given only with a 
detailed fitting of these type of sources. 


6 SUMMARY 

Following the identification of eight BL Lac objects as likely 
sources of IceCube neutrinos by PR 14, we set out to deduce the 
neutrino emission of the ones whose redshift is known by apply¬ 
ing a leptohadronic m odel, usin g a variant of the model first dis¬ 
cussed in lPetropoulou & Mastichiadisl d2012h . albeit in a more gen¬ 
eral context. We used data from near simultaneous MW obser¬ 
vations for all sources (as far as possible), and fitted the result¬ 
ing SEDs in each individual case with steady-state photon spectra, 
which were obtained with the DMPR12 numerical code. Through 
the fitting procedure we determined the distributions of injected 
electrons and protons in what we assumed were spherical volumes 
of certain radii and magnetic field strengths, moving towards us 
with certain Lorentz factors; all the above being treated as free pa¬ 
rameters. 

Of the six BL Lacs with known redshift that we studied, one 
of them, Mrk 421, had been considered for its potential as a neu¬ 


trino source in the past, by DPM14; the present study verifies that 
the neutrino event associated with it (ID 9) falls within its pre¬ 
dicted flux. A good match between modeled and observed neutrino 
(ID 22) flux was also found for 1H 1914-194. 1 RXS 1054357.3- 
553206, on the other hand, was the only BL lac to be confidently 
excluded as a source of its associated neutrino (ID 19), with the re¬ 
maining three BL Lacs displaying some discrepancy between mod¬ 
eled and observed fluxes, but not enough to rule them out as sources 
at this stage. Depending on the shape of their SEDs, we found that 
some cases favour fits dominated by SSC while others favour fits 
dominated by photohadronic interactions, with the former having 
a lower ratio of neutrino to 0.01 - 1 TeV y-ray luminosities than 
the latter. This ratio, Y vy , which could prove important for calcula¬ 
tions of the diffuse neutrino emission in the future, was then evalu¬ 
ated for its dependence on observable and model parameters. Ulti¬ 
mately, the results of our model should be confirmed or disproven 
in the near future, with the accumulation and data analysis of ad¬ 
ditional IceCube observations. Whatever the outcome, it will be 
important for constraining the available parameter space for lepto¬ 
hadronic models of blazar emission. This illustrates the importance 
of multi-messenger high-energy astronomy. 
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